{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "%reset -f\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "from os import listdir\n",
    "from os.path import isfile, join\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def processData(fileName, varindex = 1, dataindex = 10, kind = 'mol', flgerr = 'ste'):\n",
    "    working_dir = os.getcwd()\n",
    "    fileName = str(fileName)\n",
    "    callThis = join(working_dir, fileName)\n",
    "    \n",
    "    if kind == 'mol':\n",
    "        polarization_correct_mol = 0.9\n",
    "        num_pancakes_mol = 405. * np.sqrt(8/np.pi)/(1596/2 * 1e-3)\n",
    "        correct = polarization_correct_mol * num_pancakes_mol\n",
    "    else:\n",
    "        num_pancakes_na = 662. *1e-6 * np.sqrt(8/np.pi)/(0.5 * 1596e-9)\n",
    "        polarization_correct = 0.65\n",
    "        correct = polarization_correct * num_pancakes_na\n",
    "    \n",
    "#     data = np.loadtxt(callThis, delimiter='\\t', unpack=True)\n",
    "    data = np.loadtxt(callThis, delimiter='\\t', unpack=True, usecols=(varindex,dataindex))\n",
    "    \n",
    "    varindex = 0\n",
    "    dataindex = 1\n",
    "    x = data[varindex]\n",
    "    y = data[dataindex]/correct\n",
    "    x_sorted = np.sort(x)\n",
    "\n",
    "    idx_sorted = np.argsort(x)\n",
    "    y_sorted = y[idx_sorted]\n",
    "\n",
    "    var = []\n",
    "    avg = []\n",
    "    err = []\n",
    "    try:\n",
    "        for i, element in enumerate(x_sorted):\n",
    "            if i == 0: \n",
    "                i_begin = 0\n",
    "                prev_element = element\n",
    "\n",
    "            if element != prev_element or i == len(x_sorted)-1:\n",
    "            \n",
    "                i_end = i -1 \n",
    "                if i == len(x_sorted)-1:\n",
    "                    i_end = i\n",
    "\n",
    "                arr = y_sorted[i_begin:i_end+1]\n",
    "                if x_sorted[i_begin] != x_sorted[i_end] and i != (len(x_sorted)-1):\n",
    "                    raise Exception()\n",
    "                    \n",
    "                variable = x_sorted[i_begin]\n",
    "                var.append(variable)\n",
    "                avg.append(np.mean(arr))\n",
    "                if flgerr == 'stdev':\n",
    "                    e = np.std(arr)\n",
    "                else:\n",
    "                    ## this is standard error\n",
    "                    e = np.std(arr)/np.sqrt(len(arr))      \n",
    "                \n",
    "                err.append(e)    \n",
    "                \n",
    "                i_begin = i\n",
    "            prev_element = element\n",
    "            \n",
    "    except Exception:\n",
    "        print(\" --- something went wrong ----\")\n",
    "    \n",
    "    out = np.array([np.array(var), np.array(avg), np.array(err)])\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def twobody(t, initnum, gamma2):\n",
    "    return initnum/(gamma2 * t + 1)\n",
    "\n",
    "def onebody(t, amp, gamma):\n",
    "    return amp * np.exp(-gamma*t)\n",
    "\n",
    "def fullModel(n, t, gamma1, gamma2):\n",
    "    return -gamma1 * n - gamma2 * n**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "flgerr = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "list index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-10-7546502eea38>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     11\u001b[0m \u001b[1;31m# init_num = num[0]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     12\u001b[0m \u001b[0mfileName\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'978G_molonly'\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 13\u001b[1;33m \u001b[0mout\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mprocessData\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mjoin\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mworking_dir\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfileName\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;34m'.txt'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     14\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     15\u001b[0m \u001b[0mtime\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msort\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mout\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-7-2f0429666de3>\u001b[0m in \u001b[0;36mprocessData\u001b[1;34m(fileName, varindex, dataindex, kind, flgerr)\u001b[0m\n\u001b[0;32m     14\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     15\u001b[0m \u001b[1;31m#     data = np.loadtxt(callThis, delimiter='\\t', unpack=True)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 16\u001b[1;33m     \u001b[0mdata\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mloadtxt\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcallThis\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdelimiter\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'\\t'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0munpack\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0musecols\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvarindex\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdataindex\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     17\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     18\u001b[0m     \u001b[0mvarindex\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Users\\SonHyungmok\\Anaconda2\\envs\\py27\\lib\\site-packages\\numpy\\lib\\npyio.pyc\u001b[0m in \u001b[0;36mloadtxt\u001b[1;34m(fname, dtype, comments, delimiter, converters, skiprows, usecols, unpack, ndmin, encoding, max_rows)\u001b[0m\n\u001b[0;32m   1139\u001b[0m         \u001b[1;31m# converting the data\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1140\u001b[0m         \u001b[0mX\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1141\u001b[1;33m         \u001b[1;32mfor\u001b[0m \u001b[0mx\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mread_data\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0m_loadtxt_chunksize\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1142\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mX\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1143\u001b[0m                 \u001b[0mX\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Users\\SonHyungmok\\Anaconda2\\envs\\py27\\lib\\site-packages\\numpy\\lib\\npyio.pyc\u001b[0m in \u001b[0;36mread_data\u001b[1;34m(chunk_size)\u001b[0m\n\u001b[0;32m   1059\u001b[0m                 \u001b[1;32mcontinue\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1060\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0musecols\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1061\u001b[1;33m                 \u001b[0mvals\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mvals\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0musecols\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1062\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvals\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[0mN\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1063\u001b[0m                 \u001b[0mline_num\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mi\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mskiprows\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mIndexError\u001b[0m: list index out of range"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAADYBJREFUeJzt3HGI33d9x/Hny8ROprWO5QRJou1YuhrKoO7oOoRZ0Y20fyT/FEmguEppwK0OZhE6HCr1rylDELJptolT0Fr9Qw+J5A9X6RAjudJZmpTALTpzROhZu/5TtGZ774/fT++4XHLf3v3uLt77+YDA7/v7fX6/e+fD3TO/fH/3+6WqkCRtf6/a6gEkSZvD4EtSEwZfkpow+JLUhMGXpCYMviQ1sWrwk3wuyXNJnrnC7Uny6SRzSZ5O8rbJjylJWq8hz/A/Dxy4yu13AfvGf44C/7T+sSRJk7Zq8KvqCeBnV1lyCPhCjZwC3pDkTZMaUJI0GTsn8Bi7gQtLjufH1/1k+cIkRxn9L4DXvva1f3TLLbdM4MtLUh9PPvnkT6tqai33nUTws8J1K35eQ1UdB44DTE9P1+zs7AS+vCT1keS/13rfSfyWzjywd8nxHuDiBB5XkjRBkwj+DPDe8W/r3AG8WFWXnc6RJG2tVU/pJPkycCewK8k88FHg1QBV9RngBHA3MAe8BLxvo4aVJK3dqsGvqiOr3F7AX01sIknShvCdtpLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDUxKPhJDiQ5l2QuycMr3P7mJI8neSrJ00nunvyokqT1WDX4SXYAx4C7gP3AkST7ly37O+CxqroNOAz846QHlSStz5Bn+LcDc1V1vqpeBh4FDi1bU8Drx5dvAC5ObkRJ0iQMCf5u4MKS4/nxdUt9DLg3yTxwAvjASg+U5GiS2SSzCwsLaxhXkrRWQ4KfFa6rZcdHgM9X1R7gbuCLSS577Ko6XlXTVTU9NTX1yqeVJK3ZkODPA3uXHO/h8lM29wOPAVTV94DXALsmMaAkaTKGBP80sC/JTUmuY/Si7MyyNT8G3gWQ5K2Mgu85G0m6hqwa/Kq6BDwInASeZfTbOGeSPJLk4HjZQ8ADSX4AfBm4r6qWn/aRJG2hnUMWVdUJRi/GLr3uI0sunwXePtnRJEmT5DttJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNDAp+kgNJziWZS/LwFda8J8nZJGeSfGmyY0qS1mvnaguS7ACOAX8GzAOnk8xU1dkla/YBfwu8vapeSPLGjRpYkrQ2Q57h3w7MVdX5qnoZeBQ4tGzNA8CxqnoBoKqem+yYkqT1GhL83cCFJcfz4+uWuhm4Ocl3k5xKcmClB0pyNMlsktmFhYW1TSxJWpMhwc8K19Wy453APuBO4AjwL0necNmdqo5X1XRVTU9NTb3SWSVJ6zAk+PPA3iXHe4CLK6z5RlX9sqp+CJxj9A+AJOkaMST4p4F9SW5Kch1wGJhZtubrwDsBkuxidIrn/CQHlSStz6rBr6pLwIPASeBZ4LGqOpPkkSQHx8tOAs8nOQs8Dnyoqp7fqKElSa9cqpafjt8c09PTNTs7uyVfW5J+UyV5sqqm13Jf32krSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSE4OCn+RAknNJ5pI8fJV19ySpJNOTG1GSNAmrBj/JDuAYcBewHziSZP8K664H/hr4/qSHlCSt35Bn+LcDc1V1vqpeBh4FDq2w7uPAJ4CfT3A+SdKEDAn+buDCkuP58XW/luQ2YG9VffNqD5TkaJLZJLMLCwuveFhJ0toNCX5WuK5+fWPyKuBTwEOrPVBVHa+q6aqanpqaGj6lJGndhgR/Hti75HgPcHHJ8fXArcB3kvwIuAOY8YVbSbq2DAn+aWBfkpuSXAccBmZ+dWNVvVhVu6rqxqq6ETgFHKyq2Q2ZWJK0JqsGv6ouAQ8CJ4Fngceq6kySR5Ic3OgBJUmTsXPIoqo6AZxYdt1HrrD2zvWPJUmaNN9pK0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqYlDwkxxIci7JXJKHV7j9g0nOJnk6ybeTvGXyo0qS1mPV4CfZARwD7gL2A0eS7F+27Clguqr+EPga8IlJDypJWp8hz/BvB+aq6nxVvQw8ChxauqCqHq+ql8aHp4A9kx1TkrReQ4K/G7iw5Hh+fN2V3A98a6UbkhxNMptkdmFhYfiUkqR1GxL8rHBdrbgwuReYBj650u1VdbyqpqtqempqaviUkqR12zlgzTywd8nxHuDi8kVJ3g18GHhHVf1iMuNJkiZlyDP808C+JDcluQ44DMwsXZDkNuCzwMGqem7yY0qS1mvV4FfVJeBB4CTwLPBYVZ1J8kiSg+NlnwReB3w1yX8mmbnCw0mStsiQUzpU1QngxLLrPrLk8rsnPJckacJ8p60kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNDAp+kgNJziWZS/LwCrf/VpKvjG//fpIbJz2oJGl9Vg1+kh3AMeAuYD9wJMn+ZcvuB16oqt8HPgX8/aQHlSStz5Bn+LcDc1V1vqpeBh4FDi1bcwj4t/HlrwHvSpLJjSlJWq+dA9bsBi4sOZ4H/vhKa6rqUpIXgd8Ffrp0UZKjwNHx4S+SPLOWobehXSzbq8bci0XuxSL3YtEfrPWOQ4K/0jP1WsMaquo4cBwgyWxVTQ/4+tuee7HIvVjkXixyLxYlmV3rfYec0pkH9i453gNcvNKaJDuBG4CfrXUoSdLkDQn+aWBfkpuSXAccBmaWrZkB/mJ8+R7g36vqsmf4kqSts+opnfE5+QeBk8AO4HNVdSbJI8BsVc0A/wp8Mckco2f2hwd87ePrmHu7cS8WuReL3ItF7sWiNe9FfCIuST34TltJasLgS1ITGx58P5Zh0YC9+GCSs0meTvLtJG/Zijk3w2p7sWTdPUkqybb9lbwhe5HkPePvjTNJvrTZM26WAT8jb07yeJKnxj8nd2/FnBstyeeSPHel9ypl5NPjfXo6ydsGPXBVbdgfRi/y/hfwe8B1wA+A/cvW/CXwmfHlw8BXNnKmrfozcC/eCfz2+PL7O+/FeN31wBPAKWB6q+fewu+LfcBTwO+Mj9+41XNv4V4cB94/vrwf+NFWz71Be/GnwNuAZ65w+93Atxi9B+oO4PtDHnejn+H7sQyLVt2Lqnq8ql4aH55i9J6H7WjI9wXAx4FPAD/fzOE22ZC9eAA4VlUvAFTVc5s842YZshcFvH58+QYuf0/QtlBVT3D19zIdAr5QI6eANyR502qPu9HBX+ljGXZfaU1VXQJ+9bEM282QvVjqfkb/gm9Hq+5FktuAvVX1zc0cbAsM+b64Gbg5yXeTnEpyYNOm21xD9uJjwL1J5oETwAc2Z7RrzivtCTDsoxXWY2Ify7ANDP57JrkXmAbesaETbZ2r7kWSVzH61NX7NmugLTTk+2Ino9M6dzL6X99/JLm1qv5ng2fbbEP24gjw+ar6hyR/wuj9P7dW1f9t/HjXlDV1c6Of4fuxDIuG7AVJ3g18GDhYVb/YpNk222p7cT1wK/CdJD9idI5yZpu+cDv0Z+QbVfXLqvohcI7RPwDbzZC9uB94DKCqvge8htEHq3UzqCfLbXTw/ViGRavuxfg0xmcZxX67nqeFVfaiql6sql1VdWNV3cjo9YyDVbXmD426hg35Gfk6oxf0SbKL0Sme85s65eYYshc/Bt4FkOStjIK/sKlTXhtmgPeOf1vnDuDFqvrJanfa0FM6tXEfy/AbZ+BefBJ4HfDV8evWP66qg1s29AYZuBctDNyLk8CfJzkL/C/woap6fuum3hgD9+Ih4J+T/A2jUxj3bccniEm+zOgU3q7x6xUfBV4NUFWfYfT6xd3AHPAS8L5Bj7sN90qStALfaStJTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ18f+GmWq6NWLIwgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize = (6, 4))\n",
    "fileName_list = ['1288G', '980G']\n",
    "color_list = ['blue', 'red', 'black', 'darkgreen', 'orangered']\n",
    "\n",
    "working_dir = os.getcwd()\n",
    "# callThis = join(working_dir, '978G_molonly' + '.txt')\n",
    "# out = np.loadtxt(callThis, delimiter='\\t', unpack=True)\n",
    "# time = np.sort(out[0])\n",
    "# num = out[9]\n",
    "# num = num[np.argsort(out[0])]\n",
    "# init_num = num[0]\n",
    "fileName = '978G_molonly'\n",
    "out = processData(join(working_dir, fileName + '.txt'))\n",
    "    \n",
    "time = np.sort(out[0])\n",
    "num = out[1]\n",
    "num_err = out[2]\n",
    "\n",
    "# init_num = num[0]\n",
    "# num = num/init_num\n",
    "################\n",
    "## two-body fit\n",
    "# popt2, pcov2 = curve_fit(twobody, time, num, sigma =num_err, absolute_sigma=True, p0=[num[0], 1e-4])\n",
    "popt2, pcov2 = curve_fit(twobody, time, num, p0=[num[0], 1e-4])\n",
    "print(\"two body loss rate = {} (1/sec)\".format(popt2[-1]*1000.))\n",
    "twobody_coeff = popt2[-1]/popt2[0] ## this is beta constant divided by volume\n",
    "\n",
    "# ax.scatter(time, num/init_num, c = 'k', label='1288G (NaLi only)')\n",
    "ax.scatter(time, num, c = 'k', label='978G (NaLi only)')\n",
    "ax.errorbar(time, num, yerr = num_err, ls ='', c = 'k')\n",
    "plt.plot(time, twobody(time, popt2[0], popt2[-1]), c = 'k', linestyle='--')\n",
    "\n",
    "\n",
    "###################################\n",
    "def onetwobody(t, initnum, gamma1):\n",
    "    N0 = initnum\n",
    "    gamma2 = twobody_coeff \n",
    "    sol = odeint(fullModel, N0, t, args=(gamma1, gamma2))  \n",
    "    return sol[:, 0]\n",
    "###################################\n",
    "\n",
    "i = 0\n",
    "for fileName, color in zip(fileName_list, color_list):\n",
    "#     out = np.loadtxt(callThis, delimiter='\\t', unpack=True)\n",
    "    out = processData(join(working_dir, fileName + '.txt'))\n",
    "    \n",
    "    x = out[0]\n",
    "    num = out[1]\n",
    "\n",
    "    init_num = num[0]\n",
    "#     num = num/init_num\n",
    "    \n",
    "    err = out[2]\n",
    "#         err = err/init_num\n",
    "    t = x\n",
    "    ax.scatter(t, num, c = color, label=fileName_list[i])\n",
    "    if flgerr == True:\n",
    "        ax.errorbar(t, num, err, c = color, ls='')\n",
    "    \n",
    "    #######\n",
    "    ## curve_fit to the full model with ode_int\n",
    "    popt, pcov = curve_fit(onetwobody, t, out[1], sigma = err, absolute_sigma=True, p0=[out[1][0], 1e-1])\n",
    "    popt, pcov = curve_fit(onetwobody, t, out[1], p0=[out[1][0], 1e-1])\n",
    "    \n",
    "    tt = np.linspace(0, 400, 100)\n",
    "#     plt.plot(tt, onetwobody(tt, *popt)/popt[0], c = color)\n",
    "    plt.plot(tt, onetwobody(tt, *popt), c = color)\n",
    "    print(\"for \" + str(fileName) + \", loss rate = {} (1/sec)\".format(popt[-1]*1000))\n",
    "    i += 1 \n",
    "    \n",
    "ax.set_ylabel('NaLi number',fontsize = 15)\n",
    "ax.set_xlabel('Hold time (msec)',fontsize = 15)\n",
    "# ax.set_ylim([0., 1.2])\n",
    "ax.set_xlim([-5, 410])\n",
    "# ax.set_xlim([-5, 50])\n",
    "ax.legend()\n",
    "\n",
    "plt.xticks(fontsize = 12)\n",
    "plt.yticks(fontsize = 12)\n",
    "\n",
    "# plt.savefig('lifetime comparison.pdf', dpi = 500)\n",
    "# plt.savefig('lifetime comparison.png', dpi = 500)\n",
    "plt.show()\n",
    "\n",
    "########\n",
    "## simple one-body fit\n",
    "out = processData(join(working_dir, '1288G.txt'))\n",
    "popt, pcov = curve_fit(onebody, out[0], out[1], p0 = [out[1][0], 1e-2])\n",
    "print(\"loss rate from exponential fit = {} (1/sec)\".format(1000 * popt[-1]))\n",
    "\n",
    "# final = np.vstack((x_temp, num_temp, err_temp)).T\n",
    "# np.savetxt('aggregated.txt', final, delimiter ='\\t', fmt ='%f')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C:\\Users\\SonHyungmok\\Dropbox (MIT)\\BEC3-CODE\\Hyungmok\\Jupyter Notebook\\loss rate coeff_state 8 Na_aggregated\\lifetime_comparison\\for paper\n",
      "[[0.00000000e+00 1.00000000e+02 2.00000000e+02 3.00000000e+02\n",
      "  5.00000000e+02 8.00000000e+02 1.00000000e+03 1.50000000e+03]\n",
      " [3.18385528e+01 2.87664561e+01 2.36332148e+01 2.28397644e+01\n",
      "  2.16112813e+01 1.78671182e+01 1.70288829e+01 1.11238058e+01]\n",
      " [1.03957410e+00 2.15383391e+00 1.28074694e+00 1.51637899e+00\n",
      "  1.70742238e+00 7.96992340e-01 1.61683854e+00 1.21135431e+00]]\n",
      "two body loss rate = 1.00865875059 (1/sec)\n",
      "for 976G, loss rate = 80.927035174 (1/sec)\n",
      "for 999G, loss rate = 2.14355452435 (1/sec)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAESCAYAAAAIfCk9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd4FNX6wPHvSQiQhN5LSIIE6T2goICCggUVRQVBBUVRsXP1ehW9P1uu7XrtDURRCagoFrCiAoKKSJWmBIQETOgoJQkE8v7+OLubzZJNNmQ3u5u8n+eZZ9mZ2ZkzmbBvzpxz3mNEBKWUUqooEcEugFJKqdClQUIppZRXGiSUUkp5pUFCKaWUVxoklFJKeaVBQimllFcaJJRSSnmlQUIppZRXGiSUUkp5VSXYBSirBg0aSGJiYrCLoZRSYWXZsmW7RaRhSfuFfZBITExk6dKlwS6GUkqFFWNMui/76eMmpZRSXmmQUEop5ZUGCaWUUl6FfZuEUqrs8vLy2LZtG7m5ucEuivKz6tWrExcXR1RU1Al9XoOEUopt27ZRs2ZNEhMTMcYEuzjKT0SEPXv2sG3bNlq2bHlCx9DHTUopcnNzqV+/vgaICsYYQ/369ctUQ9QgoZQC0ABRQZX1vmqQUEop5VXFDBJnnGEXpVTYeO655+jYsSMdOnTg2Wefda0fPnw4Xbt2pWvXriQmJtK1a1fANraPHj2aTp060a5dOx577DHXZ3bs2MHIkSM56aST6NGjB7179+ajjz4q8rxZWVkMGTIEgPnz52OMYfbs2a7tQ4YMYf78+cWWfcyYMXzwwQeF1mVmZnLppZeW6mfgzRlnnFHsoOEjR47Qr18/jh496pfzuauYQUIpFVbWrFnD5MmTWbJkCatWrWLOnDmkpaUB8N5777Fy5UpWrlzJsGHDuOSSSwCYOXMmhw8fZvXq1SxbtozXXnuNLVu2ICIMHTqUfv368ccff7Bs2TLeffddtm3bVuS5//e//3H99de73sfFxZGSklLma2rWrNlxgSNQqlatysCBA3nvvff8fmwNEkqpUktNTSUxMZGIiAgSExNJTU0t0/HWr1/PqaeeSkxMDFWqVKF///7H/eUvIrz//vtcccUVgH3WfujQIY4ePUpOTg5Vq1alVq1afPfdd1StWpUbb7zR9dmEhARuvfXWIs/94Ycfcs4557jed+nShdq1azN37tzj9n344Yfp2bMnHTt2ZNy4cYiI12vasmULHTt2PG69iHD33XfTsWNHOnXq5Ppinz9/PmeccQaXXnopbdu2ZdSoUccdf8qUKdx5552u95MnT2bChAkADB06tMz3oSgVL0ikpsLixbBgASQm2vdKKb9JTU1l3LhxpKenIyKkp6czbty4Mn1BdezYke+//549e/aQnZ3N559/ztatWwvts3DhQho3bkzr1q0BuPTSS4mNjaVp06bEx8dz1113Ua9ePdauXUv37t19Ou/mzZupW7cu1apVK7T+/vvv59FHHz1u/1tuuYVffvmFNWvWkJOTw5w5c0p9rbNmzWLlypWsWrWKb775hrvvvpusrCwAVqxYwbPPPsu6dev4448/+OGHHwp9dsSIEXz66afk5eUB8Oabb3LNNdcA9mf4yy+/lLo8JalYQSI1FcaNg8OH7fv0dPteA4VSfjNx4kSys7MLrcvOzmbixIknfMx27dpxzz33cPbZZ3POOefQpUsXqlQpPIxrxowZrloEwJIlS4iMjCQzM5PNmzfz9NNP88cffxx37JtvvpkuXbrQs2fP47ZlZWXRsOHxiVD79u0L2MDkbt68eZxyyil06tSJ7777jrVr15b6WhctWsQVV1xBZGQkjRs3pn///q4v9169ehEXF0dERARdu3Zly5YthT4bGxvLgAEDmDNnDr/99ht5eXl06tQJgMjISKpWrcqBAwdKXabiVKwgMXEiePzykp1t1yul/CIjI6NU6301duxYli9fzvfff0+9evVcNQaAo0ePMmvWLIYPH+5aN336dM455xyioqJo1KgRp512GkuXLqVDhw4sX77ctd9LL73Et99+y65du447Z3R0tNcxBBMnTizUNpGbm8v48eP54IMPWL16Nddff/0JjT8o7hGVe40mMjKyyIbo6667jqlTpxaqRTgdPnyY6tWrl7pMxalYQcLbL2kZf3mVUgXi4+NLtd5XO3fuBGywmTVrVqFawzfffEPbtm2Ji4srdL7vvvsOEeHQoUMsXryYtm3bMmDAAHJzc3nllVdc+3rWfJxOPvnk4/5adxo0aBD79u1j1apVAK6A0KBBAw4ePHjCjdL9+vXjvffe49ixY+zatYvvv/+eXr16+fz5U045ha1btzJ9+vRCP6M9e/bQsGHDE06/4U25BwljzDRjTJYxZr8xZoMx5jrH+kRjjBhjDrotD5Tq4N5+Scv4y6uUKpCSkkJMTEyhdTExMWXuETRs2DDat2/PBRdcwEsvvUTdunVd2959991CX4hgHyMdPHiQjh070rNnT6655ho6d+6MMYaPP/6YBQsW0LJlS3r16sXo0aN54oknjjtnbGwsrVq1YuPGjUWWaeLEia5eUXXq1OH666+nU6dODB069LjHVzfccANxcXHExcXRu3dvr9d58cUX07lzZ7p06cKAAQN48sknadKkic8/J4DLL7+c0047rdDPaN68eZx33nmlOo5PRKRcF6ADUM3x77bAdqAHkAgIUKU0x+vRo4e4TJsmEhMjAgVLTIxdr5Tyat26daXaf9q0aZKQkCDGGElISJBpYfx/bNasWTJx4sRgF6NUzj//fPnmm28Krbv44ovlt99+K3L/ou4vsFR8+I4t9wR/IuLe0iOOpRWwp8wHHzXKvo4daxuvExIgJaVgvVLKL0aNGsWoCvL/6uKLL2bPnrJ//ZSHv/76i169etGlSxcGDhzoWn/kyBGGDh1KmzZt/H5OI8U0ogSKMeZlYAwQDawA+gENgM1AJjZwzAXuFpHdRXx+HDAOID4+vkd6uscsfM7R1iWMklRKWevXr6ddu3bBLoYKkKLurzFmmYgkl/TZoDRci8h4oCbQF5gFHAZ2Az2BBOzjp5pAkX1XRWSSiCSLSHJR3deUUkr5R9DmkxCRY8AiY8yVwE0i8jzgTE6ywxhzC5BljKklIvtLdXCtQSillF+EQhfYKtg2CU/O52Cly3N79CiMHm1HXSullCqTcg0SxphGxpgRxpgaxphIY8xg4ArgO2PMKcaYNsaYCGNMfeB5YL6I/F2qk6SlwVdfQe/eMGIEeOkDrZRSqmTlXZMQ4CZgG7AP+C9wh4h8ApwEfAkcANZg2ymu8HIc79q1s4HigQfg00+hTRv45z/hr7/8dQ1KqQApKl24t1ThAL/++iu9e/emQ4cOdOrUyTXg7eDBg9x00020atWKbt260aNHDyZPnhyUawp7vvSTDeWl0DgJT1u3ioweLWKMSP36Is89J3L4sPf9laqkSjtOIhBWr14tHTp0kEOHDkleXp4MHDhQNmzYUGifCRMmyEMPPSQiInl5edKpUydZuXKliIjs3r1bjh49KiIiw4cPl3vvvVeOHTsmIiI7d+6Uxx9/vByvJrSUZZxEKLRJBE5cHEydCsuWQdeucPvttqbx/vt2qJ1S6oSkptokyxER/ku2XFK6cPFIFf7111+7Ri4D1K9fn8jISDZt2sSSJUt49NFHiYiwX3ENGzbknnvuKXshK6GKHSScunWDuXPhyy8hNhaGD4dTT9VeUEqdAGey5fR0+7eWv5Itl5Qu3DNV+IYNGzDGMHjwYLp3786TTz4JwNq1a+nSpYsrQKiyqTw/RWNg8GBYsQLefBMyM+HMM+G88+DXX4NdOqXCRqCSLZeULtwzVfjRo0dZtGgRqampLFq0iI8++ohvv/32uOOmpKTQtWtXmjVrVrYCVlKVJ0g4RUbCmDGwYQM89ZTtKtu1K1x9tfeeUKmpUL26DTQ6kZGq5AKZbNlbuvCiUoXHxcXRv39/GjRoQExMDOeddx7Lly+nffv2rFq1ivz8fMAm6Vu5ciX795duuJWyKl+QcIqOhrvugk2bbO+nmTPh5JPhtttgx46C/XQiI6UKCWSyZW/pwotKFT548GB+/fVXsrOzOXr0KAsWLKB9+/YkJSWRnJzM/fffz7FjxwCb5lu0HfKEVMggccYZZ3CGM39TSerWhccft91mx4yBl1+GVq1sF9q//9aJjJTykJICHpnCiYmx68vKW7rwolKF161blwkTJtCzZ0+6du1K9+7dOf/88wF4/fXX2bNnD0lJSfTo0YOzzjqryFThqmRBSfDnT8nJybJ06dJC65wBYv6JNExv2GADxPvvQ716sHdv0fsZA47qrFLhrrQJ/lJT7d9JGRm2BqHJlkNb2CX4C6TU1FQWL17MggULSExMLP3k7CefDO+9Z7vNnnKK9/10IiNViY0aZZvw8vPtqwaIiqtCBYnU1FTGjRvHYUf7QXp6OuPGjSt9oADo3h0+/9zWKjy70kVH+6durZRSIa5CBYmJEyceN5dtdnY2E8vSfvDww/DWW9CoUcG6mjVtIsEiJilXSqmKpEIFiQwvffC8rffZlVfaHk/5+TB7NjRvbhu5O3SwD2cdPSiUUqqiqVBBIt5LO4G39aVmDAwZAkuXwocfQtWqNoB07AgzZmiwUEpVOBUqSKSkpBDj0TcvJiaGFH+3H0REwCWXwKpVthdUZCSMHAmdO8O772qwUEpVGBUqSIwaNYpJkyZRtWpVABo3bsxTTz3FyJEjA3PCiAi47DKb1uO99+y6K66wNQt9DKVUqRWVKnzVqlX07t2bTp06ccEFF7hGTh85coRrrrmGTp060aVLl0Jd3jVVuP9UqCABNlAkJCQAsGPHDm6++Wbq1KlDcnIymzZtAuCPP/7gp59+Yvfu3f4ZhRkRAZdfDqtX25pFlSr2MVT79vD229rArZQP1qxZw+TJk1myZAmrVq1izpw5pKWlcd111/H444+zevVqLr74Yp566ikA15f+6tWrmTt3Lv/4xz9cqTiuu+466tatS1paGitWrODLL79kr7cxT6pYFS5IADRp0oSePXvy2Wef8eyzz3LVVVfRoEED6tWrB8Cbb75Jnz59aNiwIXXr1qVnz56MHDmSQ4cOAZCVlcWePXtKH0CcNYtVq2ybRXS0nUq1TRuYPBmOHPH3pYa1Uo2MV6ElALnCvaUK//333+nXrx8AZ599Nh9++CEA69atY+DAgQA0atSIOnXqsHTpUk0V7mdVSt4l/Hz//ffFbr/xxhs59dRTSUtLcy2rVq1ytWfce++9vPXWW9SpU4fWrVvTunVrOnTowH333QfA4cOHqVatmvcTONsshg6FOXPgkUdsvqdHHrF5osaOtQFEqXDkzGfm7G7uzGcGZRpV17FjRyZOnMiePXuIjo7m888/Jzk5mY4dO/Lpp59y0UUXMXPmTFf68C5duvDJJ58wYsQItm7dyrJly9i6dSvbt2/XVOF+VO5BwhgzDRgIxALbgSdF5HXHtoHAS0A88DMwRkTS/V2G5s2b07x5c6/bx40bR5cuXVwB5IcffmDFihWuIHH++eezYsUKVwBp3bo13bp144ILLih8oIgIuPBCuOAC+PprGyRuvRUefRQmTIAbb4Ratfx9eUoFVnH5zMoQJNxThdeoUcOVKvyNN97gtttu4+GHH+bCCy90tTlee+21rF+/nuTkZBISEujTp0+h1OJOKSkpzJw5k507d5KZmXnC5ausyj13kzGmA7BRRA4bY9oC84HzgXRgE3AdMBt4BOgrIqcWd7yicjcFwtGjR12/gFOmTOGXX35xBZGtW7dy1llnMXfuXABOO+008vPzXQEkKSmJbt260bZNG/j+e/jPf2zQqFPHZp299VZo0CDg1xBqypRjS/lVqXI3RUQUPbOjn/OZ3XfffcTFxTF+/HjXug0bNnDllVeyZMmS4/bv06cPr7/+OlWrVmXQoEFs3LixUG2iRo0aHDx40G/lCydlyd0U1PmpgTZAFnA5MA740W1bLJADtC3uGMXOcV1OsrOzJSsry/X+lltukTPPPFPi4uIEEECuueYaERHJz8+XAQMGyP+df76sb99eBORo9epyePx4kYyMYF1CUPTv31/69+8f7GIoKeUc1wkJIjZMFF4SEspcjh07doiISHp6urRp00b27t3rWnfs2DG56qqrZMqUKSIicujQITl48KCIiHz99dfSt29f13Euu+wyuffee11zXufk5EhMTEyZyxeuyjLHdVDaJIwxLwNjgGhgBfA5kAKscu4jIoeMMZuADsBvHp8fhw0q/hsoVwbR0dFEu7UxvPDCC65/5+TksGnTJlcV+eDBgxw7dozJK1bwUGYm7YB/5uZy1auvwqRJHL78ch7KzqZGz56uWkjr1q2pUaNGeV+WUkVLSSncJgF+yxU+bNgw9uzZQ1RUlCtV+HPPPcdLL70EwCWXXMI111wD2LknBg8eTEREBM2bN+edd95xHef111/n7rvvJikpiXr16hEdHa2pwk9Q0FKFG2Migd7AGcATwKvALhH5l9s+PwCTRWSqt+OU1+OmQDh06BCbNm0iLS2NjjVr0mbOHPInTyYiN5dPsD+Unxz7Tp06ldGjR5ORkUFqaqoreCQlJYVtANHHTaGjtKnCNVd4eCnL46ag9W4SkWPAImPMlcBNwEHAsxW3FnCgvMtWXmJjY+ncuTOdO3e2KwYNIuKBB+DFF7nwxRe5aO9edrdpw7xevej4wgvw5pusnDDB1YDu1LRpU2bNmsWpp57Kxo0b+fXXX0lKSiIpKem4EehK+cWoURoUKolQ6AJbBWgFrAVGO1caY2Ld1lceDRvCQw9h/vlPmDKFBk8/zWXvvGMbBUW4cPNm9k+ezMYePVwN5xs3bnT11po9ezYTJkxwHa558+a0bt2a6dOn07RpUzZv3szBgwdJSkoq9IhMKaWK5EvDhb8WoBEwAqgBRAKDgUPARUBD4G9gGFAd+7RlcUnHDIWG64CaOlWkatXCDYRRUSKvvFLk7gcOHJBly5bJu+++K4888oiMHj1a+vTpI9nZ2SIiMmHCBFdjelxcnJx55ply/fXXS15enoiI7Nq1S3JycgJ+WdOmTZNq1aoJIAkJCTJt2rSAn1N5t27dOsnPzw92MVQA5Ofnl6nhulzbJIwxDYEPgC7Y0d7pwPMiMtmx/SzgRSCBgnESW4o7Zji3SfgkMdEOVvJkDNx0E9x5JyQl+Xy4tLQ0li1bVqgW8tdff7Fu3ToAhg8fzsyZM2nRooWrC2/nzp256aabAMjPzy/zICXn5FDuc3/ExMQwadIkRukjjKDYvHkzNWvWpH79+hhjgl0c5Sciwp49ezhw4AAtW7YstM3XNokKOcd1heKtTzrYVOV5eXDRRTZY9O1rg0cZfPXVVyxevNgVQNLS0khISGD58uUA9O3bl61btxYaA9KjRw/69+/v8zkSExNJLyLwJSQksGXLljKVX52YvLw8tm3bRm5ubrCLovysevXqxMXFERUVVWi9BomKwltNIiEBFi+Gl16Cl1+GvXvtlKt33mmTDTq63PpDdna2qwH8f//7H8uWLXMFkH379nHhhRfyySefAHD66acTGxtbqPtu586dadGihet4ERERRebFMsa4ErQppQIrLAbT+WOp8G0S06aJxMQUbpOIibHrnQ4dEnn1VZG2be32pk1FHn1UZOfOgBdv9+7dkp6eLiJ2sNPIkSMlOTlZateu7Wr7GD9+vIiI5OXlybnnnis1a9Z0bXNfEvwwGEsp5RtCsU0iECp8TQJsn/SxY+HwYVuD8NYnPT8fvvoKnnkG5s6FatXsfrffbidEKkciwu7du0lLS6Nu3bq0a9eOnTt3cu6557Ju3brjHmvExMTw+OOPM2/ePBISEkhISCA+Pp6EhATatGkTtmNBlApVlfpxk4j9voyMDFKhQsG6dfD883Y+i5wcOOMMmyfqggvsfBdBJCK88sor3HHHHeTl5dG0aVOeeuopOnfuzIgRI0hPT3elbQfb0D1y5EiWL1/OvffeWyiAxMfH0717d2rWrBnEK1Iq/FTaILF5s53rZ9IkuOqqIBYsVOzdC1OmwIsvFoyOHT8errsO6tcPatG8jbgWEfbu3UtGRgbp6en07NmT5s2bs3DhQiZMmEBGRgY7d+507b9w4UJOP/10Zs+ezRNPPHFcEOnfvz+xsbHleGVKhb6QH3EdKA0bQm4uaEZgh3r14O67bYP27Nnwwgvwr3/Bgw/aeblvuQW6dQt2KQsxxlC/fn3q169PN7ey9e3bl19++QWwObEyMjLIyMhwjVg3xlClShV++ukn3n//fY46ZgTcsmULsbGxvPTSS7z++uuFgkhCQkKh9NNKqcIqXJCoUcNO0aBBwkOVKnDxxXZZvdrWLKZNgzfegD594Oab4dJL/dorKpCio6Np06YNbdq0ca0bMmQIQ4YMAeDYsWNs376d9PR012j0Bg0a0Lx5czZt2sR3333HgQMHMMa42kfuuusuPvnkk0JtIi1btuTqq68GbA1HxxCoyqbCPW4CaNcOOnSADz4IUqHCxV9/wdSpthvtxo3QqBFcf73N8FkO2XWDmeBPRPj777/5888/6dChA2CTKH7xxRekp6eTkZFBVlYWTZs2dU1Uc8kll/Dzzz8XqoV06NDBFURyc3OpXr16uV+LUifCb20SxpjqwKfAf0Rkvn+K5z9FBYmzzrJZjH/8MUiFCjf5+bY31Msv2+lWAYYMsSO6Bw2yA/oCINSzwB4+fJjdu3e7aiKTJ0/mp59+crWVZGRk0LVrV37++WcAevTowaZNmwq1h/Tu3Zsrr7wSgN27d1OvXj2dVlOFBL82XBtj9gGXisi3/iicPxUVJJo0sX8k6+DRE5Ceblv9X38ddu6Ek06yNYtrrrE1DT8K9SBRkvz8fA4cOEDt2rUBeOWVV1i3bp0rgKSnpzNw4EA+cFRpGzVqxF9//UWLFi1cgeTcc89l+PDhAGzatInmzZtrbUSVC38HibeA/SJyqz8K50+eQSI1FUaPhmPH7BOT//xHMxqfkCNH4KOP4JVXYMECiIqCYcPghhugf/8yp/+A8A8SvsjLyyMqKgoR4bXXXiM9Pd21ZGRkMHz4cP773/+Sm5vrysrbuHFjVxAZNWoUQ4cOJS8vjzVr1pCQkEDdunW1bUSVmb+DxEjgKewcOJ8DO7CjZF1E5PMTK2rZuAeJ1NSiJ8yaNEkDRZmsXw+vvQZvvWWraCefbH/Qo0dXyrm5/cnZGJ6bm8v777/vqoE4l9tuu42bb76ZTZs2keRI5FijRg1XELnzzjs5++yz2b9/P7/++isJCQk0a9aMyEo9SEj5wt9BoqSEOiIiQfmtdA8SxaU50rxxfpCdDTNn2qj744+2J9SwYbax+4wz/FK7UEXbv38/c+fOLdQekp6ezoMPPsgFF1zAvHnzGDBgAACRkZHExcWRkJDAE088wamnnkpmZiarV692tZXoZFTK30EioaR9RKSIr+fAcw8S3hKmGmPbZpUfrV5tg8W0abZ2kZRkB+iNGQONGwe7dJXOvn37WLJkSaFHWenp6Tz33HN069aNt99+m9GjXXN60aBBA+Lj45k+fTpt2rRh/fr1rF+/3lVDadCggT7SquAq5YhrrUkEQU6O7Ws8eTIsXGjHYwwZYgPG4MFBTwGirH379rF27dpCASQjI4M33niDJk2a8J///IeJEye69o+OjiY+Pp5FixbRoEEDFi5cyIYNG2jWrJlrqV+/vvbUCmN+DxLGmGrAtUAy0AK4WUTSjDHDgV9FZH1ZCnyiSmqTiIiw6Yu0TaIc/PabHZz31lu2Z1SzZrZmce210KpVsEunirF//342bdpUqE0kIyODd999l8jISMaPH88rr7xS6DPR0dEcPHiQiIgIXnvtNVavXl0oiMTFxdG+ffsgXZEqib8fN50MzAVqA8uAM4CeIrLcGPMiUEtErvbhONWAl4GzgHrARuA+EfnCGJMIbMZOZ+r0hIg8Utwxi+rdNHGirVFERNg21vVBCV+VWF6eHW/x+uvw5Zf2WV///jZYXHqp7U2gwsrhw4fZvn07mZmZrmX//v2u2sfNN9/MjBkz2Ldvn+szLVq0ICMjA4DrrruO9evX06xZM5o2bUqzZs1o06YNF198MWDnLImOjtZHXOXI30HiSyAWuAA4CBwBkh1B4jLsl/lJPhwnFrgbmApkAOcBM4BOjl02A1EicrTEQjl4SxXepQv8+Sd07AgVuIdl6Nu2zVbl3ngDNm2CmjVhxAhbw+jdWxu7K5icnByysrLIzMwkJyeHs88+G4AHHniAH3/8kczMTLKysvj777/p3bs3PzpGvHbp0oW0tLRCNZHTTz+dW265BYAlS5ZQu3ZtmjVrphl//cTfQeIQcJmIfG6MiQTyKAgS/YCvRCT6BAv6K/AQtobityBx3nn2EXmPHhokQoKIvSFvvGF7SGVn22remDE2XW9cXLBLqMrRoUOHOHDgAE2aNAHgtddeIy0trVBNpU+fPkydOhWA+vXrs3fvXgBiY2Np1qwZV155Jf/+978BeP7552nUqFGhmopm/i2ev7PA5gLegkBz4C9fC+bOGNMYOBlY67Y63Rgj2Mdbd4vI7hM5drNmdjyYChHGQL9+dnnhBfjwQ3jzTbjvPvt88Kyz7LiLiy/Wx1GVQGxsbKEv8RtuuKHY/T/44INCASQrK4v6jlT32dnZ3H777cd95r777iMlJYVDhw5x4403FqqlOB93NdBxPiXyNUjMBe4zxnyDfdwEII42hluxA+xKxRgTBaQCb4nIb8aYGkBPYCVQH3jJsX1wEZ8dB4wDiPeSiM4ZJMK881bFVLOmrUGMGWMfQb39tl2uvNKm8b3sMrj6ahtQtPdM0IXCyPgzzzzT67bo6Gj27dtXKIBkZmZyyimnALB3714WLVpEZmYmR9z+cnzuuee47bbbSEtLY+jQoccFkXPPPZekpCTy8vLIz8+nWrVqAb/OUORrkLgb+AHb0DwXO9r630AHoCpwSWlOaoyJAN7Btm3cAiAiBwHnc6MdxphbgCxjTC0R2e/+eRGZBEwC+7ipqHM0a2ZftTYR4lq1gocegv/7P/s46q0VXWd6AAAgAElEQVS37OOoN9+0eVWuusoubinBlXJnjKFOnTrUqVOnyN5ULVq0YPPmza7JrJzB5OSTTwbsqPe2bduSmZnJ/PnzycrKIi8vjw8//JCkpCTmz5/PoEGDaNCgQaHHWXfffTft2rVjx44dbNmyhWbNmtGkSROioqLK+0cQUKXpAlsXmAAMBBoAe4Fvgf+JyB6fT2i7L7wBJALniUiOl/0aA9uBOiLyt7fjeWuTmD0bLrzQzqezfLmvpVMhITsbPv4Y3nkHvv7a9o7q2dPWNEaM8HuiQVW8UKhJlKf8/Hz27NlDbGwsMTExpKWl8e677x5XU/n4449JTk7mzTff5NprrwVswGrYsCHNmjVj5syZJCUlsXz5cpYsWVKoltKoUSOqBHkMUcgOpjPGvAp0Bc5y1B6c60/Btm2kAXWxXWUbiYj3eibeg8Sjj8IDD9h/JyRASoqOlQhLWVkwfbrt27xihZ24fNAgezOHDgVtnAy4yhYkSiszM5Ply5cXCiCZmZlMmTKFRo0a8cgjj7ga2J0iIiLYvn07DRs2ZMaMGcyfP79QEGnatCldu3Yt1WDF0t6ngExfaoypA3QEmgKZwFoR8bnR2pHe4wbgMLDdrU/0DUA+8B+gEbAf+1jritKUzyk11WZ/dUpPt4PsQANF2GnaFP7xD7usXWvTgKSm2lpFTIwNFKNGwdln20y1yq9SU1NZvHgxhw8fJjExkZSUFEbpf6JCnF/s3tx7772MHTu2UMN7Zmamq+E9LS2Njz76iF27drk+ExUVxeHDhwG45557mDdvXqEgkpCQwFVXXQXYnmLODMIBISIlLthg8gS20TrfbTkIPIntturTsfy99OjRQzwlJIjYJuvCS0LCcbuqcHTsmMiCBSLjxonUrWtvboMGIuPHiyxcaLerMps2bZrExMQItg1SAImJiZFp06YFu2gV0uHDhyU9PV0WL14sn332mWv9M888I4MHD5ZOnTpJ/fr1BZCWLVu6tg8ePFiioqKkWrVqUrNmTRkyZIhMnDjRtf2LL76QOXPmyM8//yybNm2S/fv3S35+vgBLxZfvf592guex3WDvA9piR0u3BSY61j/vy3ECsRQVJIwpOkgYU7qbpgKsf3+7lMXhwyKffCIyfLhIdLS90fHxIv/8p8iKFSL5+f4oaaWUkJBQKEA4lwT9ayuocnJyJCsry/V++vTp8q9//UsaN24sderUkS5dusgll1zi2t6hQ4fj7uGgQYN8DhKlmZnuERH5XxHb/gHcLyJ1S1uL8Yei2iQ00V+YcDxD9dtox4MH4ZNPYMYM+OorOHrU9ooaMcIubdv65zyVREREBEV9PxhjyNe0yiHHW5vEli1b2LFjB7t27WL37t3s2rXLORjRr20S+RQe8OZuDR4TEAVbSkrRkw+lpASvTKoc1Khh2ydGjYI9e2DWLHj3XXj4YdvNtnNnGD4cLr/cpjZXxYqPjye9iL+2vI1NUqEpMTGRxMTE49Y7514vkS/VDeBZYKaXbR8QYo+bRESmTROJiLBPH6pVE2nb1se6nCof06bZG+NsLArkc+7MTJHnnhPp06fg2WP37iJPPCHyxx+BO2+Y0zaJ8DFt2jSpVq2a63GgL/eIsrZJAOPdljuAbdjaxGPAnY7XdcBW4HZfThaIxVuQEBFJTLRX2Ldv2R99Kz+aNk0kJqZwg1FMTGADhVN6usjTT4v06lVw7uRkkSefFNm8OfDnDzMn8uWjyteJBnNfg4TXNgkfpiz1qJAEf/pST23awIYN0KsXREdror+QESqNRps32wmTZs6EX36x65KTbVqQYcN0DgwHHScR2hITE4t8LJiQkMCWYv4/+TpOwutIDRGJKMUSkrOuO1OtaGqOEOOYY8Dn9YHSsiXcfTcsWQJ//AFPPGFzRd1zj22z6N7dNmT99lv5lkupUsjw8v/G2/rSqtDZ06pWta+OMSkqVHhr+Axmg2jLlvDPf8LPP9vazNNPQ/XqcP/90K4ddOgA//43rFqlWSNVSPHWkcBfHQxKFSSMMW2MMQOMMed5Ln4pjZ85axIbNsCCBfYpR2pqUIukwP517pkOPJS6nyUkwIQJ8OOPsHUrPP88NGxoy9e1K7RubQPKTz/ZvFIV3Pz58/VRUwhLSUkhxuP/U0xMDCn++v/kS8MFdua4NcAxCo+4di7HfDlOIJbiGq7btCncNlqe7aOqBOXZu8lftm8XmTRJ5JxzRKKibNmbNRO56SaRr76yA/uUCoJA9m7ydTDdL9gxFfdh04Uf95RfRIpoiQy84hquq1cv+lGTDqoLEf4eTFee/voLPvvMjsX48ks7KKd2bRgyxOaTGjzYzpuhVDkJdoK/dsAwEfnKx/1DgrcG6/JuH1VehGNwcKpTp2DgXk4OzJ0LH31kc9SnptpnnQMH2oBxwQXgmKZTqXDja5BYAoTdMMv4+KJ7WuqAUeVX0dF28pILL7SpQH74waYH+fhj+NwxaeMpp8BFF9l92re307kqFQZ8bbgeB4wzxowyxjQzxsR4LoEs5IlKSQHPeT1CqX1UVUBVqkD//vC//9mpWVetgkcesQ3c990HHTva7rV33AHffQd5ecEusVLF8rVNog4wmWKmKZUQHEwHcMUVNn0P6ORDKsj+/NM+jpo9G7791jaY1a4N55xjH0mdcw445hhQKtD8OjOdMWYO0Bt4He8N12+dQDnLrKQg4ZyfpmdPO2ZKqZBw6JBtx5g9G+bMgZ077UC+006zjd/nn6+PpVRA+bvh+kzgehGZXrZilT/nhFE6oE6VC197bMXG2kbtoUPto6ilS22wmD3bjvi+5x47sOe882zAOPNM2/ahVDnztU1iC5Bd0k4lMcZUM8ZMMcakG2MOGGNWGGPOdds+0BjzmzEm2xgzzzHdaZk4g4Sm5lABl5oKixeXfuRmRIRNMPbww3Ye761b4bXXoFMnmDrVBon69W0N4+WXtf+2Kle+Bom7gYnGmMQynq8KNmtsf6A28ADwvjEm0RjTAJjlWFcPWAq8V8bzaU1ClY/UVDuJifMXzTmx+okM8Y+Ls5/99FM7L8YXX8B118H69XDzzTaFSPv2dt7vb77RX24VUKUZTBcP1MXWKv7y3EdEep1QAYz5FXgIqA+MEZE+jvWxwG6gm4h4zbBWUpsE2DFNY8fCs8+eSAmV8kF5ZLYVgbQ02632889tjeXIEfvoasAAOPdc2/jdsqV/zqcqtDJngfWwBvgcSAV+wM4r4bmcSCEbAyc7Pt8BWOXcJiKHgE2O9WXSvDls21bWoyhVjPLIbGsMnHyy7T779dewd6+tbVx9NaxeDePHw0kn2Wla77jD1kCyy/iUODXVBsCICE1+Vkn5VJMIyImNiQK+ADaJyA3GmCnALhH5l9s+PwCTRWSqx2fHYcduEB8f36OoXOruzj0XduyAWrXs+3Ae6KtCVLDnyBCxmSy//NIu8+dDbq4d+d2vn00TMniwzWbra48p5yM0z3mAJ03SfuQVgL9rEn5ljIkA3sF2pb3FsfogUMtj11rAAc/Pi8gkEUkWkeSGDRuWeL6kJFtL1wzPKmCCndnWGDvL1u232xrE3r3w1Ve2DSMzE+66yzaEN28OY8bA9Om2221xJk48viaSnW3Xq0rDpy6wxpj3S9pHRC738VgGmAI0Bs4TEeeQ07XAaLf9YoFWnOCjLHetW8PBg3Zwq3OOCaX8yvmX9dixtiE52CM3o6Nh0CC7PP207TE1d659TDV7NrzlGNbUrRucfbZdTj/dZsV0CpXJoVRQ+dpwPa+I1fWANsAe4HcRGeDTCY15FegKnCUiB93WN8QO1LsW+AzbmN1fRE4t7ni+NFx/8YXtbh4VZQNFsP//KhVUx47ZrrZffWUDx48/2v8Y1atD375w1ll2ufjiogOCplGuEPw6mE5EzvRykhbAR8AzPhYqAbgBOAxsNwXPRm8QkVRjzDDgRWAa8DMwwpfjlmTdOvvqTJPj7J0IGihUJRQZaefyTk62j44OHoTvv7cB45tv7EA+gBo17L7HjhV8VpOfVTplbrh2fLE/KiLt/FOk0vGlJpGQoH8QKeWzrCybW8qZNmTfPrs+MtI+kho3zo4Ab9o0uOVUZeLX3E0lnGgo8I6IBGWGFV+CRERE0Y3WxlSK2SeVOnEi8PvvNmh8+y3Mm2cnXAI7oG/AALv07w/16gW3rKpU/J3gr30Rq6tiJyN6BMjwtU3C33wJEsHunahUhXHsGKxcadOcf/stLFxoezwZA1262IBx5pm2baN27WCXVhXD30EiHyhqRwP8AlwhIn+UupR+4EuQSE21vf6OHi1Yp929lfKDI0dseuV58+zy44+2d1dEBPToYRMennmmfUyl07mGFH8Hif5FrM4FtonInydQPr/xJUiAHZT6zjv239q7SakAycmxSQ7nz7dBY/Fi22MkMtIGjf797XL66VrTCLJya5MINl+DhLMbbNeutvefUqocZGfb2sWCBXZxBo2ICPuf0Rk0+vbVNo1y5u/5JJwHrQY0B6p7bhORdaU5VnlLSrKvOTnBLYdSlUpMTMG4C7BBw5lOfcECm/r8GUcP+k6dbAoRZ9Bo0iR45VYuvo64bg68Bpxb1GZse0VQpi/1VWKirfEOGxbskihVicXEFPSIAptfaskS2wD+/fd2/oyXXrLbWre2wcK5nHSS/2bq83VyKOVzTWIy0B2YAKyjiOlLQ11UlA0UGzcGuyRKKZfq1W3toV8/O7AvLw+WL7dBY+FC+PhjeOMNu2/TpjZYnH66XTp3tn/5qYDyNUichp2+tMQcTqGsdWub6E8pFaKiouCUU+xy1112INP69QVBY9EieN/xNVSzJvTubecFP/10+5nY2JLP4ZxB8PBh+5ej9mIplq9BYicQ9k/zk5Lghx/s+CCdX16pMBARYdObd+gAN95o12Vk2IDxww92efBB+586MtI2hp92WsHSvHnh43mbQRA0UHjhaxfYK4CbsVlb9we8VKXga+8mgOeft5mUt2+Hxo0DXDClVPn46y/46aeCoPHzzwU9VOLjbbDo08cumrTQxd+9my7BTl+a7pjK1HP6UhGR4aUsY7lr3dq+pqVpkFCqwqhTx84sdq6jX01enh0V/sMPNnh8/z3MmFH8MTT9uVe+TjrUADuV6EogCmjosTQKSOn8zNkNVhuvlarAoqKgZ087het779m5i9PTbaDwNuo7JsZ2x12xonBqBlW2VOHhxtkNVhuvlapk4uPtcuzY8VOyRkZClSp2Fj+wASM5GU491S6nnALNmgWn3CGgVIPpwp12g1WqkvM2g+DIkba28dNPtk1j8WI7yM85CU1cXEGvq1NOsSlGfOlJVQFUmrQcTueeCzt22K7YSqlKypfBdLm5tm3j558Llj8ceUwjIqBjRxswevWyS/v2tkYSJgKSlqMi0G6wSimfVK9e8MjJadcuO0J8yRIbND74ACZPtttiYmwNo2fPgsWfo8SDpNyDhDHmFmAM0AmYISJjHOsTgc3AIbfdnxCRR/x5/tat4cAB2LlTezgpVWmdaDqOhg3h/PPtAvavzU2bCgLHkiW2ATw3126vV68gYDinjPUcuxHiglGTyAQeBQYD0UVsryMiAete4OwGu3GjBgmlVBkZYx9PJCXZdg2w7Rhr1sAvv9hlyRJ47LGCucKbNi0IGD162NcQ/jIq9yAhIrMAjDHJQFx5n9/ZDXbsWPjtt/I+u1KqwouKgm7d7OIczZ2dDatWFQSOZctgzpyCeZWbN7cBw7mEUODwGiSMMe8D94rIJse/i+PPwXTpxhgB5gJ3i8huPx0XsL2bQFOGK6XKUUyMzTPVu3fBugMHbMP40qV2WbYMZs8uCBzNmhUEje7d7WvTpuXexlFcTaIhduAc2MFyge4GtRvoiR2wVx94CUjFPpYqxBgzDhgHEB8fX6qTREXZ9igNEkqpoKpZsyANutOBA3ZA37Jldlm+vHCNo1EjGzC6d7c1le7doWXLgAaOoHWBNcY8CsQ5G66L2N4EyAJqF5cvqrRdYMG2JeXl2fuhlFIh7eBBW+NYscIGjeXLYd26gpHhtWvbxIbOR1zdukHbtvYv4mKUWxdYY0w/4EERGVDWY3lwRi+/h8joaNi/X7vBKqXCQI0aBXNoOOXm2sbx5ctt8FixAl59taBXVbVqdhxH164FS+fOUKtWqU/vj4brhkB/X3c2xlRxnDcSiDTGVAeOAj2wiQPTgLrA88B8EfnbD2UsJDbWdjT44w9o1crfR1dKqQCrXr2gh5TT0aOwYYMNGM6ax8cfw5QpBfu0alUQNHwUjC6w9wP/5/b+SuAh4HfgP9j2j/3Yhusr/H3y1NSCjMBt29reTuvX+/ssSilVzqpUsaO+27cvSD8iAn/+aYPGypW2h9XKlfDhhz4ftsxtEsaYYcD7IhKUeQRL0ybhnG/EPbdXRAS8/bbON6KUqkQOHMDUquVTm4SvqcIrhIkTCwcIsLMjTpwYnPIopVRQeEuZXoTixkmM9/EYvj/cCjJv84rofCNKKVW04tokXizFccIilWx8vM0GXNR6pZRSx/P6uElEIkqxBKU9orRSUuzAR3fG2PVKKaWOV6lShXvONwJQt642WiullDelarg2xlQxxpxkjGnvuQSqgP42apRND9+2ra1F7N1r8zmlpga7ZEopFXp8qkkYY6Kwg9tGA9W87BYWj5zAzky3YUNBOpT09IJkjVqrUEqpAr7WJP4NDAHGYtNk3AJcA3wLbAEuCEThAmXzZtv11V12tnaFVUopT74GicuBBwFnyvAlIvK2iAwCFgEXBaBsAeNsj/CkXWGVUqowX4NEC2CDiBwDcrG5lZxSgWH+LlggVfPywEy7wiqlVGG+BoksoI7j35uBfm7bwi5FXsuWNh2Hu5gY7QqrlFKefA0S8wHnzBiTgfuMMdONMW8CTwOfBKBsAbN+vc3X5KxRRETAa69po7VSSnnydZzERKABgIg8a4wxwKVANPAC8HBgihc4o0bB5Mk2QeLGjdCvX8mfUUqpysanICEi24Htbu+fAZ4JVKHKkzPP1bJl2iahlFKeKlUW2KLExkJkpA0SSimlCisuC+x3pTiOiMhAP5Sn3EVGQocOGiSUUqooxT1u2uPD55sCfQiTLLDedO8On32mc14rpZQnr0FCRC7zts0YEw/cgx2FvZtStE8YY24BxgCdgBkiMsZt20DgJSAe+BkYIyJFJPf2r0WLYNcu2LYNWrQI9NmUUip8lDbBX5IxZgqQBlwI3AskiMhjpThMJvAo8IbHsRsAs4AHgHrAUuC90pTvROzYUTDHRPfumuhPKaXc+RQkjDEdjDHTgfXAmcDtQCsReVZEckpzQhGZJSIfc/zjrEuAtSIyU0RysWlAuhhj2pbm+KVx/fU2FUdenn2/e7dN9KeBQimlrGKDhDGmhzFmFvAr0A24DmgtIq+KyBE/l6UDsMr5RkQOAZsc6z3LNc4Ys9QYs3TXrl0nfMKi5rzWRH9KKVXAa5AwxnwBLAFaAiNEpJ2IvOXI3xQINYC/Pdb9DRw3Y7eITBKRZBFJbtiw4QmfUOe8Vkqp4hXXu2mw47UF8JIx5qXiDiQijcpYloNALY91tYADZTyuVzrntVJKFa+4IPFQuZXCWoud1AgAY0wsNnng2kCdMCXFtkG4P3KKiNBEf0op5VRcF9iABAljTBXHeSOBSGNMdeAo8BHwlDFmGPAZdqKjX0Xkt0CUA46f87pmTfs6LKwSnyulVOAEIy3H/UAO8C/gSse/7xeRXdh5KVKAfcApwIhAF2bUKMjNtQPpZsyAI0dg4cJAn1UppcJDuQcJEXlQRIzH8qBj2zci0lZEokXkDBHZUp5lO+MMqFoVvvqqPM+qlFKhq9In+HMXGwunn65BQimlnDRIeBg8GNasgczMYJdEKaWCT4OEh0GD7OvXXwe3HEopFQo0SHjo3BmiouCf/wx2SZRSKvg0SHiIiIB69WDfPjgWqLHlSikVJjRIFKFuXTh6FFasCHZJlFIquDRIeEhNhT/+sP8+6yzNCKuUqtw0SLhJTbVpOo448tv+/bdNJ66BQilVWWmQcFNU6vCcHE0drpSqvDRIuNHU4UopVZgGCTfeUoTHxZVvOZRSKlRokHCTkgIxMcevP+ec8i+LUkqFguLmk6h0PFOHx8fb+a/Xrw9uuZRSKli0JuFh1Cg49VTo39/OWnfHHbBoEfz+e7BLppRS5U+DRAmuvhoiI+HNN4NdEqWUKn8aJIowf75dAJo0gfPPh7fesqOwlVKqMgm5IGGMmW+MyTXGHHQsQX/Qc+21sH07fPllsEuilFLlK+SChMMtIlLDsbQJdmHOOw8aNYLXXw92SZRSqnyFapAIKVFRNj3Hp5/CqlXBLo1SSpWfUA0SjxljdhtjfjDGnBHswgDcdRfUqQP33RfskiilVPkJxSBxD3AS0ByYBMw2xrRy38EYM84Ys9QYs3TXrl3lUqg6deBf/4LPP4fvvy+XUyqlVNAZEQl2GYpljPkS+ExEXihqe3JysixdurRcypKTA0lJdpDdjz+CMeVyWqWU8jtjzDIRSS5pv1CsSXgSICS+jqOj4cEHYfFi2z6hlFIVXUgFCWNMHWPMYGNMdWNMFWPMKKAf8FWwy+Z0zTVw8sm2bUKnN1VKVXQhFSSAKOBRYBewG7gVGCoiQR8r4VSlik0EuG4dvP12sEujlFKBFfJtEiUpzzYJJxHo1Qt27IDffis6c6xSSoWyitQmEXKMgf/+F7ZutQkAlVKqotIgcYL697ddYidPhnffDXZplFIqMDRIlMHDD0OfPjBuHGzcGOzSKKWU/2mQKIOoKJgxwzZmjxhhJypSSqmKRINEGcXH27kmli2De+4JdmmUUsq/NEj4wUUXwW23wXPPwfTpwS6NUkr5j85x7SdPPgkrVsBVV0F+Plx5ZbBLpJRSZac1CT+pVg2++ML2err6anjjjWCXSCmlyk6DhB/FxsJnn8GgQTB2LLz6arBLpJRSZaNBws+io+Hjj2HIELjpJnj8cfv4SSmlwpEGiQCoXh0+/BAuvxzuvRcGD4Y//wx2qZRSqvQ0SARI1ap2JPakSXbuiU6d4P33g10qpZQqHQ0SAWSMnRt75UqbXnz4cBg1CtLTg10ypZTyjQaJctC6NSxaZCcs+uAD+/6mm2yCQKWUCmUaJMpJlSrwf/9nczyNHQtTptipUG++GdavD3bplFKqaBokylmLFvDKK5CWBmPG2DaL9u0hORmefx527gx2CZVSqoAGiSBJSIDXXoNt2+CZZ2w32dtvh2bNYMAAO/vdzz/D0aPBLqlSqjILuZnpjDH1gCnAIOwUpveKiNeMSMGYmS5Q1qyB1FQ7cnvVKruuVi04/XTo3h26doVu3aBlS9sorpRSJ8rXmelCMUjMwNZwxgJdgc+APiKytqj9K1KQcLdrF8ybB99+a7vQrl8Px47ZbTVr2vaMVq0Klrg4aNrU1kQaNIAIrSMqpYoRlkHCGBML7AM6isgGx7p3gD9F5F9FfaaiBglPOTmwdq1NIvjrr7Bpk102b4a8vML7VqliA0X9+lCvnl3q1rW1kpo1C5aYGLtER9vX6tVtDqqqVQteo6IKL1WqQGSkfdVApFT48jVIhFoW2JOBY84A4bAK6B+k8oSM6GjbuJ3scUuPHbPtGn/+CVlZdsnMtDWRvXvtsnmzne/iwAG7+DNNSGSkXSIiCl7dF2OOfy1qgeP/7curO891vjySC5XHdqFSDqU8hVqQqAH87bHub6Cm+wpjzDhgHEB8fHz5lCxERUbaRvCEBN/2F7G1kgMHIDvbLjk59jU3186ud+RIwWteXuHl2DHbmO7+mp9f+FXE/tu5ON+7v7ovznK5/9uXV8/rKu69t59FKAiVcqjKZd063/YLtSBxEKjlsa4WcMB9hYhMAiaBfdxUPkWrGIwpeMyklKq8fK29htpT5Q1AFWNMa7d1XYAiG62VUkoFVkgFCRE5BMwCHjbGxBpjTgMuAt4JbsmUUqpyCqkg4TAeiAZ2AjOAm7x1f1VKKRVYodYmgYjsBYYGuxxKKaVCsyahlFIqRGiQUEop5ZUGCaWUUl5pkFBKKeVVSOVuOhHGmF1AUROCNsBmka1I9JrCg15TeKjs15QgIg1L2insg4Q3xpilviSvCid6TeFBryk86DX5Rh83KaWU8kqDhFJKKa8qcpCYFOwCBIBeU3jQawoPek0+qLBtEkoppcquItcklFJKlZEGCaWUUl5VuCBhjKlnjPnIGHPIGJNujBkZ7DKVljFmvjEm1xhz0LH87rZtpOO6DhljPjbG1AtmWb0xxtxijFlqjDlsjJnqsW2gMeY3Y0y2MWaeMSbBbVs1Y8wbxpj9xpjtxpgJ5V54L7xdkzEm0RgjbvfroDHmAbftoXxN1YwxUxy/UweMMSuMMee6bQ+7e1XcNYX5vZpmjMlylG2DMeY6t22Bu08iUqEWbHrx97BToZ6Onf60Q7DLVcprmA9cV8T6DthZ+vo5rm868G6wy+vlGi7BZvN9BZjqtr6B455cBlQHngIWu21/DFgI1AXaAduBc4J9PSVcUyIgQBUvnwvla4oFHnRcQwQwxPE7lhiu96qEawrne9UBqOb4d1tH2XoE+j4F/cID8MtxBDjZbd07wOPBLlspr8NbkPgPMN3tfSvH9dYMdpmLuZZHPb5QxwE/etyzHKCt4/2fwCC37Y8QYoGwiGsq6Ysn5K/Jo7y/AsMqwr0q4poqxL0C2gBZwOWBvk8V7XHTycAxEdngtm4VNgKHm8eMMbuNMT8YY85wrOuAvR4ARGQTjqAYhPKdKM9rOARsAjoYY+oCzdy3E173L90Ys80Y86YxpgFAuF2TMaYx9vdpLRXkXnlck1NY3itjzMvGmGzgN2yQ+JwA36eKFiRqYKtd7v4GagahLGVxD3AS0Bzb73m2MaYVFeP6iruGGm7vPbeFst1ATyABW/2vCaQ6toXNNfJUjecAAAhjSURBVBljorDlfktEfqMC3Ksirims75WIjMeWpy92qufDBPg+hdzMdGV0EKjlsa4W9nlk2BCRn93evmWMuQI4j4pxfcVdw0G397ke20KWiBwEljre7jDG3AJkGWNqESbXZIyJwD6aPQLc4lgd1veqqGuqCPdKRI4Bi4wxVwI3EeD7VNFqEhuAKsaY1m7rulC4mhmOBDDY6+jiXGmMOQmohr3ucOF5DbHYtpW1IrIPW4Xu4rZ/ON4/5whVEw7XZIwxwBSgMTBMRPIcm8L2XhVzTZ7C6l55qILjfhDI+xTsBpgANOi8i+3hFAucRpj1bgLqAIOxvRSqAKOAQ9iGqg7AfmxVMxaYRgg2qjmuo4rjGh7D/jXnvJ6GjnsyzLHuCQr3xHgcWIDtidHW8QseKr1LvF3TKY77EwHUx/aumxcO1+Qo36vAYqCGx/pwvlferiks7xXQCBiBfXwU6fiOOARcFOj7FPSbGYAfZj3gY8cPMAMYGewylbL8DYFfsNXBvxy/6Ge7bR/puK5DwCdAvWCX2ct1PIj9K819edCx7Sxsw1sOtidXotvnqgFvYIPhDmBCsK+lpGsCrgA2O+5JFvA20CRMrinBcR252EcTzmVUuN6r4q4pXO+V43thgeM7YT+wGrjebXvA7pPmblJKKeVVRWuTUEop5UcaJJRSSnmlQUIppZRXGiSUUkp5pUFCKaWUVxoklFJKeaVBQimllFcaJFRAGWMeNMbs9rJtqjFmaVHbijmec9KYISXsd4sxpthBQMaYQcaYO/xRrkAwxsQ4JonpH+yy+MIYE22M2WmM6Rvssij/0SChKrNBwHFBAptvf0z5FqVItwKbRWRBsAviCxHJAV7A/vxUBaFBQikPIrJJRNYEswyODKY3Y9MphJOpQD9jTKdgF0T5hwYJFVKMMV2NMd865urdZ4xJdUwaU9xnqhljXjTG/GWM2WuMeQaIKuEzDwL/ABIcj6/EOOat9nzcZIwZ49je3dj5x7ONMSsd72MdE9f8bYz5w5HW3fNcFxk7N3au4/HRk455DoozADufyCyPY4kx5k5jzNPGmD2Oianucmwb7SjDX445jau7fa6OMeZ1Y0ymoxwZxpjJHsfuaIz5zNh5oQ8YY2YaY5p47FPfGPOasXMt5xpjfnd/ZCciW7G5x64u4fpUmKho80moEGWMKep3zXjs0xCbnGw9NpFhDWwGy7nGmGQROeLl8I8D1wETgXXA9dj5fovzOtAa+2V8sWPdrhI+8xbwIjbL5uPAB8AS7CxglwLXAm8bYxaKyDbHNV2OzUr8GnAfNoXzY9g/0O4q5lwDgQ0isqeIbf8APsMmqxsCPGWMaYSdTOc2IB54BptC/nHHZ/4H9AHuxM5x3AI7VzqOciYBP2DnWrgKm2n0EeyEV71ERIwx0dj70wh4CJtQLsmxuPsRm3BOVQTBzm6oS8VeKDpzqvuy1G3fx7FZLmu5revl2O8Kx/tEx/shjvf1sZkv73H7TAT2C0xKKNt/gS1FrJ/qUa4xjnOOdlt3nmPdG27ragN5wE2O9wZIB970OP61jjLXL6ZsXwMzi1gvFE5tHYHNZrrP4+f2PvCz2/s1wK3FnO8d4Hegqtu61sAx4HzH+xuAfKBrCT/XMcBRoHqwf/90Kfuij5tUefgb+1eu5zLHY79ewNcist+5QkSWAFuA070cuxM2h/4nbp/Jd3/vR9+6/Xuj4/U7t/P+ja2NNHesOhn7V/37xpgqzsXxmepAx2LO1QQ71Wax5XBc62ZgmfvPzVG+5m7vVwJ3G2PGG2OKmhP9LOAjIN+tnJuxP/tkxz4DgBUisrKYcuModyQ2vbUKcxokVHk4KiJLPRfA81FKU2y+e087sPOEFMX5zHynx3rP9/7wl9u/jxSxzrne2RbQwPH6ObaG4Vw2O9a3KOZc1bHzF5dUDuc5iysH2Ok7Pwb+DfxujEkzxoxw294AO7d6nsdykls562NrLSVxlrt6sXupsKBtEiqUZGGfd3tqzP+3d/egTUZRGMf/x81BcRAFBb8GcSoOCoKTiIiTCrq7VKgFF7E4tm46uLlUKIJI8QspIn5AQbBSP6hCK6Ii6KAiaKsgFgfxcThvbKh522CDSeX5LQnJvXlPspzce87LhZGSOR+Kx2XARNXrtT7nX6vEcxB4UuP91zVeq567pFGBSPpC1isOR0Qb0AWcj4hRSc+K610lazXTVVY04/xZf6ilEvfEjKNsXnCSsFbyAOiIiEWSvgJExGayDjFUMmeMPIFsN1mHqLSP7q7jetP/bTfaC+AdeUrYmdkG15i7tvEhgaTRiDhKntS2gSz2D5LbXyOSym5CHAT2R0SbpNEZLrEGGFftorvNM04S1kpOAR3ArYg4wVR30xhwpdYESeMR0Qv0RMQP8oD39mLubJ4DyyPiAFnY/STpzVy/RFVsPyPiCHAuIhYDN8jEtA7YA+yTNFky/R6wNyIWFHWHOYmIIXKl8JQsfreTR3g+LIZ0F8+vR0QfuXpYCewAzkq6Qx712QncLlqIK4lsvaRjVZfbRHY42X/ANQlrGZI+AtvIlUE/cBq4S57xXdb+Crl10kfut/cD78mEM5uLZCfTSbK3v/svQy8l6QK5qtkIXCLvezgEPGaqrlHLALAQ2NqgUIbJrqPL5PdeCuxS0aor6SWwBZgEesmE1kPWF14VY76TxetrwPFiTBf5ewO/W523U5LUbf7xGddmLSoiBoC3kjqbHUu9ImInmYRWSPrW7Hhs7pwkzFpUUY8ZBFZL+tzseOoRETeB+5K6mx2LNYa3m8xalKRH5HbOqmbHUo/ijuxh6tvqs3nCKwkzMyvllYSZmZVykjAzs1JOEmZmVspJwszMSjlJmJlZqV9CovGlqmsCLAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize = (6, 4))\n",
    "# fileName_list = ['1288G', '980G']\n",
    "# fileName_list = ['976G', '981G', '999G', '1008G']\n",
    "fileName_list = ['976G', '999G']\n",
    "color_list = ['blue', 'red', 'darkgreen', 'orangered', 'black']\n",
    "\n",
    "m = os.getcwd()\n",
    "working_dir = join(m, 'for paper')\n",
    "print(working_dir)\n",
    "\n",
    "# callThis = join(working_dir, '978G_molonly' + '.txt')\n",
    "# out = np.loadtxt(callThis, delimiter='\\t', unpack=True)\n",
    "# time = np.sort(out[0])\n",
    "# num = out[9]\n",
    "# num = num[np.argsort(out[0])]\n",
    "# init_num = num[0]\n",
    "fileName = '978G_molonly'\n",
    "out = processData(join(working_dir, fileName + '.txt'))\n",
    "    \n",
    "print(out)\n",
    "time = out[0]\n",
    "num = out[1]\n",
    "num_err = out[2]\n",
    "\n",
    "# init_num = num[0]\n",
    "# num = num/init_num\n",
    "################\n",
    "## two-body fit\n",
    "popt2, pcov2 = curve_fit(twobody, time, num, sigma =num_err, absolute_sigma=True, p0=[num[0], 1e-4])\n",
    "# popt2, pcov2 = curve_fit(twobody, time, num, p0=[num[0], 1e-4])\n",
    "print(\"two body loss rate = {} (1/sec)\".format(popt2[-1]*1000.))\n",
    "twobody_coeff = popt2[-1]/popt2[0] ## this is beta constant divided by volume\n",
    "\n",
    "# ax.scatter(time, num/init_num, c = 'k', label='1288G (NaLi only)')\n",
    "ax.scatter(time, num, c = 'k', label='978G (NaLi only)')\n",
    "ax.errorbar(time, num, yerr = num_err, ls ='', c = 'k')\n",
    "plt.plot(time, twobody(time, popt2[0], popt2[-1]), c = 'k', linestyle='--')\n",
    "\n",
    "###################################\n",
    "def onetwobody(t, initnum, gamma1):\n",
    "    N0 = initnum\n",
    "    gamma2 = twobody_coeff \n",
    "    sol = odeint(fullModel, N0, t, args=(gamma1, gamma2))  \n",
    "    return sol[:, 0]\n",
    "###################################\n",
    "flgerr = True\n",
    "i = 0\n",
    "for fileName, color in zip(fileName_list, color_list):\n",
    "#     out = np.loadtxt(callThis, delimiter='\\t', unpack=True)\n",
    "    out = processData(join(working_dir, fileName + '.txt'), varindex = 1, dataindex = 10)\n",
    "    \n",
    "    x = out[0]\n",
    "    num = out[1]\n",
    "\n",
    "    init_num = num[0]\n",
    "#     num = num/init_num\n",
    "    \n",
    "    err = out[2]\n",
    "#         err = err/init_num\n",
    "    t = x\n",
    "    ax.scatter(t, num, c = color, label=fileName_list[i])\n",
    "    if flgerr == True:\n",
    "        ax.errorbar(t, num, err, c = color, ls='')\n",
    "    \n",
    "    #######\n",
    "    ## curve_fit to the full model with ode_int\n",
    "#     popt, pcov = curve_fit(onetwobody, t, out[1], sigma = err, absolute_sigma=True, p0=[out[1][0], 1e-1])\n",
    "    popt, pcov = curve_fit(onetwobody, t, out[1], p0=[out[1][0], 1e-1])\n",
    "    \n",
    "    tt = np.linspace(0, 400, 100)\n",
    "#     plt.plot(tt, onetwobody(tt, *popt)/popt[0], c = color)\n",
    "    plt.plot(tt, onetwobody(tt, *popt), c = color)\n",
    "    print(\"for \" + str(fileName) + \", loss rate = {} (1/sec)\".format(popt[-1]*1000))\n",
    "    i += 1 \n",
    "    \n",
    "ax.set_ylabel('NaLi number',fontsize = 15)\n",
    "ax.set_xlabel('Hold time (msec)',fontsize = 15)\n",
    "# ax.set_ylim([0., 1.2])\n",
    "ax.set_xlim([-5, 310])\n",
    "# ax.set_xlim([-5, 50])\n",
    "ax.legend()\n",
    "\n",
    "plt.xticks(fontsize = 12)\n",
    "plt.yticks(fontsize = 12)\n",
    "\n",
    "# plt.savefig('lifetime comparison.pdf', dpi = 500)\n",
    "# plt.savefig('lifetime comparison.png', dpi = 500)\n",
    "plt.show()\n",
    "\n",
    "########\n",
    "## simple one-body fit\n",
    "# out = processData(join(working_dir, '1288G.txt'))\n",
    "# popt, pcov = curve_fit(onebody, out[0], out[1], p0 = [out[1][0], 1e-2])\n",
    "# print(\"loss rate from exponential fit = {} (1/sec)\".format(1000 * popt[-1]))\n",
    "\n",
    "# final = np.vstack((x_temp, num_temp, err_temp)).T\n",
    "# np.savetxt('aggregated.txt', final, delimiter ='\\t', fmt ='%f')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
